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A nonlinear aircraft model is presented and used to develop an overall unified approach 
to online trim and maneuverability envelope estimation with uncertainty quantification 
without any requirement for active input excitation. The concept of time scale separation 
makes this method suitable for the adaptive characterization of altered safe maneuvering 
limitations based on aircraft performance after impairment. The results can be used to 
provide pilot feedback and/or be combined with flight planning, trajectory generation, 
and guidance algorithms to help maintain safe aircraft operations in both nominal and 
off-nominal scenarios. 
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Set of real n- vectors 

Set of n x n symmetric matrices, Set ofnxn positive definite symmetric matrices 

Trace operator, i.e., tr A is the sum of the diagonal elements of A 

Determinant, Euclidean norm for vector inputs 

True airspeed, flight path angle, course angle 

Magnitude of net thrust force produced by engines 

Magnitude of lift, drag, and side forces 

Continuous dynamics model 

State vector, Collection of state observations 

Virtual input vector, Collection of virtual input observations 

Set of allowable or achievable states and input vectors 

Average current and next state, i.e., x(k) = [x(k + 1) + x(k)\/2 

Angle of attack, Side-slip angle, Roll angle 

Acceleration due to gravity, Mass of the aircraft, Surface area of the aircraft wings 
Air density, Dynamic pressure (q = pV 2 / 2) 

Multivariate normal distribution with mean /i and covariance matrix E. 

Probability density function for the variable x given y 

Expansion coefficients for the non-dimensional coefficient of lift, drag, and side force 
Aerodynamic coefficient parameter vector, e.g., c = [Do, D\, D 2 , To, L\, Yf[ 
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State process noise precision (a.k.a., inverse covariance) matrix 
Optimal estimates for c and S 


All units are SI. 


I. Introduction 

The safe maneuvering envelope is a fundamental property of the aircraft’s design and overall current 
state of health. By definition it determines the overall capability of the aircraft. If actively monitored in 
time, it may function as an early warning system as well as provide anticipatory guidance to help avoid loss 
of control. For example, automated planning tools may use it to help pilots land safely under emergency 
landing conditions, 1-3 or when combined with a display it may provide better pilot awareness of the aircraft 
state when an automation system switches off. Additionally, the physics-based maneuverability envelope 
can be analyzed separately from the control strategy, and knowledge of the envelope may for example 
unmask control limitations clouded by adaptive controllers, and even lower barriers to the introduction of 
more advanced unconventional control strategies. 4 For these reasons, improved methods for tracking aircraft 
maneuverability in real-time may effectively help pilots avoid inappropriate crew response and further prevent 
or recover aircraft from upset conditions. 

Generally, the maneuvering envelope is the set of safe aircraft state and control inputs. Unfortunately, 
because of the underlying nonlinear aircraft dynamics, this set is difficult to calculate accurately and rapidly 
enough to provide the pilot or automation system with reliable information in a diverse and rapidly chang- 
ing environment. Previous research has considered a wide variety of approaches to meet the challenge. The 
most straight-forward of course is through wind tunnel testing, flight testing, and high-fidelity model-based 
computation. 5 8 A more sophisticated analysis is obtained by formulating the envelope estimation prob- 
lem as a reachability problem, 9 and then solving the associated Hamilton- Jacobi equations, often through 
the use of level set methods. 10 Further extensions on this theme include leveraging the concepts of time- 
scale separation 11 or semi-Lagrangian level sets . 12,13 Other alternatives rely on linearization and region 
of attraction analysis, 14 quaternion-based control architectures, 15 robustness analysis, 16 frequency domain 
non-parametric system identification, 17 or the Control- Centric Modeling methods suggested by Boeing. 18 
There even exist methods based on an artificial immune system paradigm. 19 21 

This article develops an integrated method complementary to much of the previous work on this topic. 
We begin by developing a simplified but effective nonlinear flight dynamics model in section II with important 
benefits including the ability to numerically calculate the corresponding trim envelope in real-time. The trim 
envelope however depends on aerodynamic parameters, and in section III we introduce a reliable Bayesian 
system identification technique, with a conservative approach to uncertainty quantification, tailored to the 
task of estimating the aerodynamic parameters from measured flight data. In section IV, we show how 
the dynamics model and trim envelope are combined to determine an extended safe maneuvering envelope, 
robust to the uncertainty established by the system identification process. Finally, main applications for the 
overall approach are summarized in section V. 

II. Dynamics Model 

In this section, a focused approximate nonlinear model for the aircraft dynamics is presented. This model 
captures the slower aircraft dynamics, and ultimately enables the nonlinear analysis of maneuverability in 
flight path angle and total airspeed that is central to the discussion in this article. 

Aerodynamic forces are most clearly organized around the velocity vector expressed using the kinematic 
axes, which are defined by the basis vectors 

V = cos 7 cos yx# + cos 7 sin yi/# — siny£# 

7 = — siny cosy &e — siny sin xVe — cosy £# 

X = ~ s in X &e + cos x Ve , 
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where 7 is the flight path angle, x * s the course angle, and xe,Ve, %e are unit vectors in the earth axes 
aligned with north, east, and down, respectively (see Fig. 1). In this kinematic coordinate system, an aircraft 
moving with true airspeed V has the velocity v = VV, from which one obtains the acceleration vector 

V = VV + Vj t V = VV + V 77 + VXCOS 7 *. (1) 

With this expression, one can rather easily sum up the forces and obtain the equations of motion using 
Newton’s second law. 

In summary, the forces acting on the aircraft are gravity, thrust, and the aerodynamic forces of lift, drag, 
and side force. Further assuming that the force of gravity is aligned with ze, that thrust is entirely aligned 
with the aircraft body axis xb, and following the standard convention for the directions of the lift, drag, 
and side forces, one obtains from Newton’s second law the equations of motion: 


V = 

(. Ft cos a cos j3 — Fd — mg sin 7) /m 

(2) 

7 = 

{Fd,. cos if - F &y sin^ - mg cos 7) /mV 

(3) 

X = 

(- Fa L sing? + Fa Y cos ip) /mV cos 7. 

(4) 


In these equations, Ft is the net thrust produced by all engines, Fd is the magnitude of the drag force, a is 
the angle of attack, /3 is the sideslip angle, g is acceleration due to gravity, and m is the mass of the aircraft. 
Finally, 


Fa L — Ft sin a + Ft 
Fa Y = Ft cos a sin /3 + Fy , 

are the net force magnitudes aligned with the lift and side force direction vectors, while Ft and Fy are 
the respective magnitudes of the lift and side force components. Fig. 1 illustrates the relationships between 
direction vectors and angles according to standard convention. 



Figure 1: (a) 3D orientation of velocity vector v with respect to the earth fixed axes, (b) Roll angle 
<p rotation to get from kinematic axes to aerodynamic axes, where V is into the page, (c) Side-slip 
angle (3 and angle of attack a rotations to go from the velocity axis V to the body axis x b ■ 


II. A. Maneuverability Model 

Next, parameterized expressions for the aerodynamic forces Fd, Ft , and Fy are developed. In general such 
expressions would depend on many parameters (a, /3, mass, air density, wind speed, etc.). However, one can 
obtain a simplified, but perhaps more useful model, by focusing on the primary dependencies through the 
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model defined with 


Fd{ol) — q S Cd{ol) = q S (Dq + LL<a + L^tr 2 ) 
Fl(ol) = gS Cl (a) = gS (I/ 0 + Lia) 

M0) = q S Cy(a) = q S (Yq + ^0) , 


( 5 ) 

(6) 
( 7 ) 


where q is dynamic pressure (q = pV 2 / 2, where p is air density), S is the net wing surface area, and {Di c }j: =0 
and {L/c, Yj C }J, =0 are the expansion constants for the non-dimensional aerodynamic force coefficients Cd, Cl, 
and Cy- Also note, for an undamaged aircraft one usually sets Yq = 0. Next inserting (5)-(7) into (2)— (4) , 
and relying on the small angle approximations sin/? « 0, sin a cos /3 « 0, and cos a cos /? « 1, one obtains the 
following coupled pair of nonlinear differential equations for the state vector x = [V, y] T 


V 


Xi 


—ci kx\ — c 2 nu 2 x 2 — C3 hiu^x 2 + ui/rr\ — g sin x 2 

. i 


_ x 2 


C4 KXi COS Us + C5 COS Us — CQ KU4X1 SU1 Us — g cos X2 / X\ 


(8) 


where u = [Ft, cq </?, /?] T is the virtual input vector, 

c = [D 0 , L>i, F 2 , L 0 , Li, Yi] t 

is the aerodynamic parameter vector, and where n = Sp/2m. Furthermore, the remaining possibly time vary- 
ing parameters (g, m, S, p) are assumed effectively know from measurement or separate estimation process. a 
We refer to such a model as a maneuverability model because of its dependence on virtual inputs rather than 
direct control inputs, which if needed can be determined by leveraging the principle of time-scale separa- 
tion. 11 The virtual inputs are however the parameters most closely related to any overall determination of 
maneuverability for the slower aircraft dynamics covering the motion of the center of mass. Furthermore, 
this model has a variety of properties extremely useful to applications, which we will bring to light in much 
of the discussion to follow. 

II. B. Trim 

The set of trimmable aircraft states is defined by {x \ f{pc,u\c) = 0, (x,u) E S}, where B represents the set 
of overall allowable states and virtual inputs. In addition, B should limit (x, u) to the domain over which 
the model is valid. A common approach for discovering trimmable states is to solve some variation of 

u = argmin||/(x,u;c)||, (9) 

u£t3 

for any desired fixed x, using standard constrained optimization routines. 22,23 While this approach can 
work for almost arbitrarily sophisticated aircraft models, its application to map the entire trim envelope is 
numerically intensive, and requires special care when the local minima is not equal to zero. However, in 
stark contrast to this local optimization based approach, the specific maneuverability model (8) enables a 
far better analytical approach, capable of mapping the (V, y)-trim envelope in tens of milliseconds 15 for any 
fixed aerodynamic coefficient vector c. This is an important capability because the set of trimmable states 
represents an a-priori safe aircraft maneuverability envelope. 

Characterizing the set of trimmable points for fixed c, then, requires setting the top and bottom equations 
on the right hand side of (8) equal to zero. The bottom equation is solved for u<i (angle of attack) in terms 
of the other variables in that equation, which do not include u\ (thrust). The top equation is then solved 
for u\ and we substitute the previous solution for U 2 . The result is a closed form expression for the required 
thrust and angle of attack needed to achieve trim for any given state and remaining virtual inputs. This 
enables a fast numerical sweep to determine the typically non- convex trim envelope as follows: 

a In our implementation, we consider ( g , m,S,p) as part of the virtual input vector u. 
b Using current commercially available laptop computational capability. 
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1. Setup a grid of state values in B. For most practical applications a coarse resolution is sufficient. 

2. Fix values for the roll angle p = 113, and the side-slip angle (3 = 114. 

3. For each point in the state grid, solve for the thrust u\ and angle of attack U 2 needed to achieve trim. 

4. Return only those trim points for which Ft and a are within B. 

A crude visualization of the (V,7)-trim envelope is then obtained by simply plotting the trimmable points 
from the above calculation. In addition, one should check whether the achieved trim points are stable. 
This requires checking the eigenvalues associated with the local linear approximation to (8) at each trim 
point [24, Thm. 2.1 pg. 36] and as before a closed form expression can be found, but it is too unwieldy to 
include here. The entire computation is fast enough to enable dynamic re-computation of the trim envelope 
as aircraft conditions change, or to compute extended envelopes by sweeping over values for p and /3. 

As an example, we set up a scenario using the Research Civil Aircraft Model (RCAM) 25 in landing 
configuration (flaps and gear down, air density at sea level, etc.). Table 1 shows the parameter values used 
for the setup. The actual trim envelope calculation for a high resolution grid covering 409149 points was 
completed in 114 ms on a 2.6 GHz MacBook Pro. The calculation included a check for stability, and in this 
case, all of the trimmable points were stable (but in general this may not be true). The calculated envelope 
is shown in Fig. 2. 

Table 1: Parameter values for an example maneuverability envelope calculation. 


Constants 

Aero. Coefs. 

Bounds 

Grid Res. 

S = 260 m 2 

m = 120 x 10 3 kg 

g = 9.81 m/s 2 
p = 1.225 kg/m 3 

D 0 = 0.1599 
£>1 = 0.5035 

D 2 = 2.1175 
I/O = 1-0656 

Lx = 6.0723 
U = -1 

V e [50, 150] m/s 
7 e [-20,20] deg. 

F t e [20546, 410920] N 
a e [0, 14.5] deg. 

Vres = 0.2 m/s 
7res = 0.05 deg. 


Required Thrust [kN] 
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Figure 2: Maneuverability trim envelope for a representative commercial aircraft (similar to a B767 
or A300) in landing configuration with p = /3 = 0. The plot on the left shows the maneuverability 
shaded according to the thrust required for trim, while the plot on the right shows the same envelope 
shaded according to required a. 
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In Fig. 2, it can be seen that the trim envelope boundaries are imposed by the input saturation limits. 
The upper boundary corresponds to maximum thrust Fr max = 410920 N, and the lower boundary is imposed 
by minimum thrust Fr min = 20546 N. The range for the angle of attack causes the left (<a m ax = 14.5°, prior 
to stall) and right (<a m in = 0°) boundaries. It can also be observed that more thrust is needed for larger 
flight path angles, since the thrust force has to counteract the aircraft weight, and that a larger angle of 
attack occurs for slower speeds. Analyzing the range of thrust values over airspeed for constant flight path 
angle, shows that more thrust is needed for a further speed decrease below 69 m/s. This region corresponds 
to the range of the angle of attack a > 4.5°. Analysis of the Lift-Drag Polar and the power required curve, 
shown in Fig. 3, confirms that this region is the back side of the power curve, and that V = 69 m/s is the 
minimum drag airspeed. 


Lift Drag Polar 



CD [-] 

(a) Lift-Drag Polar 



V [m/s] 



V [m/s] 

(b) Power required curve 


Figure 3: The Lift-Drag polar and power requirement curves of the RCAM model confirm that the 
minimum drag airspeed is V = 69 m/s, and the region to the left is the back side of the power curve. 


II. C. Damage Scenarios 

Damage can be injected into our maneuverability model through two primary mechanisms. First, damage 
leading to diminished control authority is modeled through a constriction of the set B of overall allowable 
states and virtual inputs. For example, the maximum permissible angle of attack may diminish due to 
icing, or the maximum available thrust may be reduced due to engine damage. The second mechanism is 
through changes to the aerodynamic parameter vector c, and this will likely occur in conjunction with the 
first mechanism. This occurrence may for example model a decrease in the conversion between a and the 
generation of lift with perhaps an additional cost in drag. 

To illustrate, Fig. 4 shows the effect of reduced aerodynamic capability, modeled by a 20% decrease to lift 
and 20% increase to drag, along with a 50% reduction to available thrust. For these damage scenarios, it can 
be seen in the left plot of Fig. 4 that a 20% decrease in lift combined with a 20% increase in drag results in a 
shift of the trim envelope towards higher airspeeds and lower flight path angles. This is physically explained 
by the force equilibrium equations: 


C L (a)pV 2 S/2 = Wcos^kW (10) 

F t -C D (a)pV 2 S/2 = Wsin'yfssWr (11) 

Due to the reduced lift capability, (10) shows that a higher airspeed is needed to compensate for the aircraft 
weight W = mg. On the other hand, an increase in drag means that reduced net excess thrust is available 
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Figure 4: Trim envelope calculation for nominal and damaged cases with (p = f3 = 0. The damaged 
case on the left corresponds to a 20% decrease of all the lift coefficients {£&}, and a 20% increase of 
the drag coefficients {Dk}. The plot on the right shows the same damage case with an additional 50% 
loss of available thrust. 


for climb in (11). If the maximum thrust limit Fr max is reduced by 50%, then there is a further reduction 
in net excess thrust available for climb, which is translated into a shift downwards of the upper limit of the 
trim envelope, as shown on the plot to the right in Fig. 4. Furthermore, this 50% reduction to FV max makes 
establishing trimmed level flight impossible in this scenario. 

The maneuverability model we just developed may not accurately capture the aerodynamics for many 
other damage scenarios. In this instance, alternate parameterizations of the aerodynamic force equations 
(5)-(7) may be required, and existing research is in place to support this task. 26 28 These extended cases 
can be treated by a system identification procedure capable of using observed flight data to select between 
competing damage models, like the one developed in the next section. 

III. System Identification 

The maneuverability model (8) is characterized by a vector of aerodynamic parameters c, which once 
known, ultimately determines the trim envelope as discussed above. A critical feature of the same ma- 
neuverability model, is that it is linear in these parameters for given state and virtual inputs. Therefore, 
the application of almost any system identification algorithm, say from [29], should yield quite satisfactory 
results with sufficient input excitation to cover the full dynamic response of the aircraft. The question of 
input or maneuver design for system identification is important and covered in the Refs. [26-31], along with 
additional methods for developing more sophisticated aerodynamics models and online parameter retrieval. 

However, many important applications for online maneuverability characterization are aligned with track- 
ing nominal performance, and identifying degradation to the envelope caused by impairment, during normal 
operations that would not constantly explore the space of possible inputs. Therefore, in an effort to enable 
progressive learning of the aircraft dynamics with the possibility of insufficient excitation, a novel Bayesian 
method was developed with the following high-level features: 

• works in the background from pilot inputs and measurements collected while maneuvering naturally, 

• leverages optimization and probability theory to simultaneously estimate both system parameters and 
uncertainty from the measured data, 

• is biased towards greater uncertainty quantification, especially when excitation is insufficient, 

• seamlessly enables selection between alternate impairment models through the evidence calculation. 
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The objectives of this section, then, are to develop fast methods for estimating the aerodynamic coefficient 
vector c from measured state and input data, to quantify the resulting estimation uncertainty, and to enable 
the comparison between competing models that explain the data. 

III. A. Measurement Model 

Our system identification algorithm stems from a midpoint-Euler discrete approximation for the measured 
state process: 

x(k + 1) = x{k) + hf(x , u ; c) + r(fc), (12) 

where h is the time-sample resolution, x = [x(k) + x(k + l)]/2 and r(k) ~ AT(0, hiS -1 ), for k = 0, 1, . . . , m — 
1. While the maneuverability model is linear in c, it is still nonlinear in x and u, and the midpoint- 
Euler approximation enables better representation of the true nonlinear system, for negligible additional 
computational cost. Furthermore, the more general case is developed where x(k),r(k) G M n and c G and 
this makes the method readily applicable to broader more sophisticated dynamics models. Clearly, n — 2 
and d = 6 for the specific maneuverability model discussed above. Finally, the inverse covariance matrix S 
is an important statistical parameter that characterizes the process noise r(fc), which we will also estimate 
from the observed data. This matrix is sometimes also referred to as the precision or information matrix, 
and we use these terms interchangeably in the discussion to follow. 

III.A.l. Sensor Measurement Error 

A more comprehensive approach would also model the effect of additive measurement noise to the states 
and inputs: 


y x (k) = x(k)+r] x (k) 

y u {k) = u(k)+r] u (k), 

where rj X:U may be an independent additive Gaussian process. However, we will avoid this consideration by 
assuming the variance of r(k) far exceeds that of rj XjU , which is not an impractical assumption in many cases. 
With additional sophistication, our approach is extendable to include the measurement noise as well. 


III.B. Probabilistic Approach 

Our objectives are met by using the measurement model (12) to develop an expression for the posterior 
probability density function (pdf) 

p(c,5|x,u,x(0)), (13) 

where x is short-hand for [x(l),x(2 ), . . . , x(m)\ and u is short-hand for [i/(0), 'u(l), . . . , u(m — 1)]. This is the 
pdf that describes the joint set of reasonable aerodynamic parameters c and process noise precision matrices 
S', given the available data (x, u, and x(0)). The remainder of this section outlines the basic derivation of 
the posterior pdf. Then, starting in §III.C, a fast algorithm to find the estimates for c and S that maximize 
this posterior distribution is presented, along with analytical approximations for quantifying the estimation 
uncertainty. 

Using Bayes’ Theorem the posterior distribution is decomposed as follows: 

p(c S\x u ic(0)) = p( c > S,x|u,3;(0)) = p(x|c, S, u, a;(0 ))p(c, S|u, s(0)) ^ 

where we refer to p(x|c, S, u, x(0)) as the likelihood distribution, p(c, Sju, x(0)) as the prior distribution, and 


Z = f p(c, 5, x|u, x(0)) dc'dS' 

Jc',S' 


(15) 


as the evidence. The evidence Z = p(x|u, x(0)) is the probability density of the observed state data given the 
input measurements, initial state, and all implied modeling assumptions. In essence, it is the probability of 
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the observations given the model. This is an important quantity because the evidence for alternative models 
can be computed and compared; a process that would allow one to extend the present framework to consider 
alternate maneuverability models that may arise in damage scenarios, as previously discussed in II. C. 

III.B.l. Priors 

The prior distribution is assumed to take the following form 

p(c,5|u,x(0)) =p(c,S) = p(S)p(c). (16) 

This distribution characterizes the available information prior to making any measurements, and the equation 
implies prior independence of the parameters, inputs, and initial state. 

The prior on the aerodynamic parameter vector c is further assumed to follow a multivariate normal 
distribution, 

p(c) = (27r) _d/2 |S c |“ 1 / 2 exp|-i(c-/Lt c ) T E“ 1 (c-/i c )|, (17) 

where fi c and £ c are used to set the prior knowledge of the model parameters. 

The prior on the precision matrix S is assumed to follow a specialized version of the Wishart distribution: 

p{S) = ^„ + iexp|-itr(A _1 S')| , (18) 

where T n+ i is the known normalizing constant of integration, and A E §++ is a given matrix that enables an 
encoding of prior information as needed. 0 Typically, one would set A to a diagonal matrix related to a prior 
worst case assessment of expected model error and turbulence. This creates a pessimistic prior on S that 
forces a justification of uncertainty through measurement. To see how this works, consider the case where 
the airspeed and flight path angle errors are independent and expected to be Ai and A 2 , respectively. Note, 
these are error variances based on how well one expects the model to perform before any new measurements 
are taken; and if this is unknown one should set Ai and A 2 to cover the range of possible values (usually 
known from the underlying physics of the problem). Setting A = diag(Ai, A 2 ), the Wishart prior looks like, 

p(S)oce x p(-|l)exp(-|i), (19) 

where S\% and S 22 are the diagonal elements of S. This, however, is equivalent to the case where Sn and 
S 22 are each independently distributed according to an exponential distribution with means 2Ai and 2 A 2 - d 
Furthermore, exponential priors assign higher probability density to small values for Sn or S 22 , and since 
S is an inverse covariance matrix this corresponds to larger variances in the expected model performance. 
This is what we mean by pessimistic prior, and it is a key feature of our approach. 

The ability to leverage prior information in the probabilistic formulation enables a number of distinct 
advantages especially when used in combination with evidence calculations. For example, in a nominal 
aircraft one may have a high degree of confidence in the values of the aerodynamic parameters c from flight 
and wind tunnel testing, which can be encoded into p(c). In this case, the observed data will support the 
model with a sharp prior and the corresponding evidence calculation will be high. However, in the presence 
of damage the aerodynamic parameters will shift, contradict the prior, diminish the evidence calculation, 
and force the consideration of alternative models. 

c This corresponds to the general Wishart distribution with n + 1 degrees of freedom. 

d In this case, there is also no prior preference on S12 other than the implicit constraint required for S to be positive definite 
(he., S'nS '22 — Sf 2 > 0 ==> —VS11S22 < S21 < VS11S22). 
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III.B.2. State Likelihood 


Next we develop the likelihood function for the states given the aerodynamic parameters, precision matrix, 
virtual inputs, and initial state. Following Ref. [32], this is accomplished by first noting the system of 
equations induced by (12), represents a nonlinear transformation of random variables, from r(k) to x(k + 1) 
for k = 0, . . . , m — 1. The corrected likelihood pdf is thus given by: 


m-i , 1 

^(x|c, S, u, x(0)) = |J| (27r/i) _n/2 |S'| 1/2 exp] ~—b(k) T Sb{k) 

k = 0 ^ 


( 20 ) 


where, 


b(k) = x(k + 1) — x(k) — hf(x(k),u(k ); c), 


( 21 ) 


and where J is the mn x mn Jacobian matrix accounting for the full transformation. For small h, the 
determinant of J is well approximated by 


| J\ ~ exp 


, m — 1 

2 Y ^F(x(j),u(j);c) 


where F(x,u;c) is the n x n Jacobian matrix of the vector field / with respect to x. One should notice that 
tr F = V • /, which is the divergence of the vector field defined by /. 

Many nonlinear dynamical models, including the maneuverability model (8), can be expressed as affine 
functions of the unknown model parameters [32, §II.D] e , generally taking the form 


f(x, u; c) = U(x, u)c + v(x, u), (22) 

and when this is the case one can also write, 

tr F(x,u]c) = q(x,u) T c-\- r(x,u). (23) 


Thus, the nonlinear dynamics are characterized as a combination of the columns of U(x,u), with unknown 
weights c, along with a known part specified by v(x,u). 


III.C. Estimation Algorithm 


With all the required distributions now defined we have two remaining tasks. The first is to find the maximum 
a posteriori (MAP) estimates for the aerodynamic parameter vector c and precision matrix S by solving 

maximize lnp(c, S^x, u, x(0)) = lnp(x|c, 5, u, x(0)) + lnp(c) + In p(S) — In Z , 


where lni^ is a constant that can be ignored for this step. The second task is to determine an approximation 
for the marginal posterior distribution p(c|x, u, ^(o))- 

Finding the MAP estimate is equivalent to minimizing the negative log of the joint pdf p(c, S', x|u, ^(0)) 
with x fixed by solving: 

minimize /o(c, S) = — In p(c, 5, x|u, *(0)) (24) 

with respect to the variables c eR d and S G §+ + and where, 


m— 1 


fo(c, S) = ^ In 2t rh - ^ In |S| + | ^ <gc + r k + ( U k c + v k ) T S ( U k c + v k ) 


k=0 


(25) 


+ ^ln \(2-rr) d T, c \ + i(c- ^ C ) T S C 1 (c - fi c ) + t tr(A 1 S) - ln'I'n+j. 


In this expression, = q(x(k),u(k)), rk = r(x(fc), iz(fc)), Uk = U(x(k),u(k)), Vk = v(x(k),u(k)) — Ax(k), 
and A x(k) = [x(k + 1) — x(k)\/h. One way to accomplish the minimization (24) is by using the following 
block coordinate descent algorithm: 


e If an exact representation is not possible, then usually the nonlinear model can at least be approximated in this manner. 
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1. Start from an initial guess (c°, S' 0 ), and k = 0. 

2. Precision Update. S /c+1 = argmin 5e § n + fo(c k ,S). 

3. Parameter Update. c /c+1 = argmin c /o(c, S /c+1 ). 

4. Repeat until Ac T MAc < e, where Ac = c k+1 — c fe , and M is defined by (26) below. 

It can be shown that each update step requires solving a well-posed convex optimization problem and 
this condition is sufficient to guarantee convergence of the overall algorithm to a joint local stationary 
solution [33, §2.7, pg. 272]. Furthermore, by separately setting the gradient with respect to S and c equal 
to zero, one finds, respectively, that the precision update has a closed form analytical expression, and 
the parameter update requires solving a well posed regularized least-squares problem. While we can not 
guarantee global joint optimality in both c and S', the algorithm does uniquely determine the global optimal 
precision matrix corresponding to any fixed estimate for c (and vice-versa), and in practice the algorithm 
has demonstrated rapid convergence to reliable estimates. 

In keeping with the probabilistic approach, one may desire to stop the algorithm when Ac is insignificant 
with respect to the uncertainty in the current parameter estimate c k . This prevents needlessly updating the 
current estimate when the uncertainty is large compared to the update magnitude. One can show that 

M=(x- 1 +hJ2 Ul S k+1 U k J , (26) 

V k=0 J 

specifies the approximate multivariate precision matrix for the current parameter estimate. Therefore, the 
stopping criterion that Ac T MAc < e, ensures the update falls within a small uncertainty ellipsoid around 
the current parameter estimate. In fact, one can show that Prob {Ac t MAc < e} = F x 2 (e), where F x 2 is the 
cumulative distribution function (cdf) of a chi-squared distribution with d degrees of freedom. In particular, 
we used e = 0.001, for which F x 2(e)«2 X 10 11 , in our examples/ 

Finally, an important consequence of the convexity inherent to both the precision and parameter updates, 
is that a variety of hard constraints can be incorporated into the problem. In particular, if the constraints 
are also convex then the overall convergence properties of the algorithm are unaffected (see Refs. [33] and 
[34] for further details), although in some cases more careful analysis of the estimation uncertainty may be 
required. For example, known trim points are easily handled by introducing linear equality constraints into 
the parameter update step. 


III.D. Parameter Uncertainty 


Determining the estimation uncertainty amounts to finding a quadratic approximation for the negative log 
of the joint distribution and normalizing appropriately. This is readily accomplished by forming a quadratic 
second order taylor approximation around the optimal estimate. For (25), the general form is 


fo(c,s) = f 0 (c*,S*) + ± 



Vg/o V cs /o 

V cs /o T V 2 Jo 



(27) 


where the gradient term is omitted because it is approximately zero at the MAP parameter estimate (c*, s*). 
Furthermore, S is expressed through the expansion 

z 

sT s i £ i- ( 28 ) 

3 = 1 

f To find the value of e that produces, say, a 0.1% confidence ellipsoid we simply lookup or numerically compute e = 
F~^ (0.001). For example, this is accomplished in Matlab® using chi2inv(0 . 001 , d). 

x d 
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where the Ej G § n serve as basis matrices for the space of possible precision matrices. 8 After substituting, 
(28) into (25), the expressions for V^/o, V cs /o, and V^/o can be found. In particular, one may readily 
observe that V^/o = M, where M is defined by (26) evaluated with 5' /c+1 = S'*. 

With these expressions, we find from (27) that the approximate full posterior distribution is 


p(w |x, u, ^(0)) 


|(27r) d+2 S|- 1 / 2 exp 


{ 



where w = (c — c*,s — s*) G R d+Z , and with H = —V 2 In |5|, 


V 2 / 0 V cs /o 

V c ,/o T f H 


(29) 


(30) 


The approximate marginal posterior distribution for the aerodynamic parameter vector is then, 

p(c|x,u,a:(0)) = |(27r) d £'r 1/2 exp |-^(c - c*) T S ,_1 (c - c*)| , (31) 

where E /_1 = V;?/o — (2/m)V cs /oi7 _1 V cs /q\ Asymptotically, as m oc, we see E /_1 « V^/o = M, and 
observe in (26) that M is essentially the prior information matrix updated by the data information matrix. 
Also, for large m we have that p(c |x, u, x(0)) ~ _p(c|*S''*', x, u, x(0)). That is, the marginal posterior for the 
aerodynamic coefficient vector is approximately the conditional distribution evaluated at the precision matrix 
estimate. Furthermore, p(c|5*, x, u, x(0)) = p(c|S*, x, u, x(0)) exactly, because /o is quadratic in c. 


III.E. Evidence 

The evidence calculation is derived from (27) as follows. The approximation for the joint pdf is 
p(c,5,x|u,a:(0)) = exp {-/ 0 (c, S)} 

~ exp(-/o(c,s)| 


= exp{-/ 0 (c*,S'*)}exp|-iu; T E 

= p(c*,S*,x u..r(()))exp ( -I(r 7 


The evidence is then obtained by integrating the approximate expression over all possible values for w 
The result in log-form is 


In Z w lnp(c*,5*,x|u,a:(0)) + bn|(27r) <i+z E| 
= — /o(c*, S*) + ^ In |(27r) d+2 S|, 


(c, s). 


where /o is given by (25). The accuracy of this evidence estimate can be checked using non-deterministic 
methods such as importance sampling [35, §13.3]. 


III.F. Extensions 

The system identification method developed so far is amenable to a variety of extensions. An important 
example is the incorporation of direct accelerometer measurements. This can be accomplished, for example, 
by introducing an accelerometer measurement process model defined by, 

a(k) = A k c + r a , (32) 

g While all of § n is covered by z = n{n + l)/2 such basis matrices, there are implicit (convex) constraints on the feasible 
values s can take, since S must also be positive definite. 
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where a(k) = [Fd, Fl, Fy] T /m represents the component aerodynamic accelerations (corrected for gravity 
and thrust) observed at each time instance &, A k is the corresponding matrix from (5)-(7) defined as 


A k = KXx (k) 2 


1 a(k) a(k) 2 0 0 0 

0 0 0 1 a(k) 0 

0 0 0 0 0 /3(k) 


(33) 


and r a is the associated measurement noise vector distributed Af(0 , S'” 1 ). Here, S a is a separate unknown 
precision matrix for the accelerometer error. Recall from before that x\ (k) is the total airspeed (which is 
the first element of the state vector), and k = Sp/2m. Now, the full joint distribution can be factored as 


p(c, S, S a , a, x|u, x(0)) = p( a, x|c, S, S a , u, x(0)) p(c, S, S a |u, x(0)) 

= p(a|c, S a , X, u, x(0)) p(x|c, S, u, x(0)) p{c) p{S ) p(S a ), 


where a represents the full set of accelerometer measurements, analogous to x and u. The second line in 
this factorization follows by applying the definition of conditional probability to factor p( a, x|-), and because 
the prior knowledge regarding c, S, and S a is assumed to be assembled from independent sources before 
accounting for any new measurements. A similar approach to the one outlined in §III.B.1-III.E can be 
followed to derive the estimation algorithm, this time with an additional Wishart prior characterized by A a 
(compare (18)) and update step for S a , along with the calculations for uncertainty and evidence. 11 

Finally, we wish to emphasize the overall formalism developed here is readily applied to dynamical sys- 
tems that are nonlinear in the state and input variables, but affine in the unknown system identification 
parameters. This characteristic leads to the analytical uncertainty and evidence approximations just dis- 
cussed, which in turn make it possible to quickly evaluate competing models using standard techniques from 
Bayesian model selection [36, Ch. 4]. We will highlight one such example in the discussion to follow. 


III.G. Example 

The overall system identification approach just discussed is now demonstrated using an example that focuses 
on verifying our estimation algorithm and uncertainty quantification under the condition that all of our 
modeling assumptions hold. To accomplish this, we started by simulating the full nonlinear dynamics for 
a general transport class aircraft, similar to a B757, operating nominally at an altitude of 3048 meters 
for 90 s. This was done only to obtain realistic values for the virtual inputs. Next, using these inputs, 
we generated state measurements exactly from our model (8) in time step increments of 0.01 s, while also 
adding independent zero- mean Gaussian state noise with standard deviations of 3 m/s and 0.5° for V and 
7, respectively. Accelerometer measurements were also generated concurrently with the state measurements 
according to the model (32), with accelerometer measurement covariance S = 0.1 2 /. The prior distribution 
on the parameter vector was set using 

Me = [0, 0.5, 0, 0, 4.5, 0] T £ c = diag (2 2 , 3 2 , 3 2 , 3 2 , 3 2 , 2 2 ) , 

where diag(-) denotes a diagonal matrix with diagonal elements defined by the input list. Also, the Wishart 
priors for S and S a were characterized by diagonal matrices A and A a , based on assuming worst case 
standard deviations of 10 m/s for V, 5° for 7, and 1 m/s 2 for the accelerometer measurements. Finally, 
to simulate a rapid fault onset, the first 45 s of our simulation corresponded to the nominal model with 
true underlying coefficients from Table 1, while the final 45 s were produced using coefficients from the 
20% damage model discussed in §11. C. Fig. 5 shows the angular virtual inputs, state, and accelerometer 
measurements corresponding to our simulation. 

Our estimation procedure was then applied to the above simulated data after decimation by a factor of 
10, so that the measurements correspond to a 0.1 second time sampling resolution. This decimation step was 

h The particular simplistic approach suggested here assumes independent state and accelerometer measurements. This can 
be realized, for example, by using GPS data to derive the state measurements separately from the accelerometer measurements. 
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Figure 5: Left Column: virtual angular and state inputs to the estimation procedure produced using 
a general transport aircraft simulator. Right Column: measured accelerometer data and model fit 
at the estimated parameter values. The onset of a fault occurs 45 s into the simulation. The fault 
corresponds to 20% increase in all drag coefficients and a 20% decrease in all lift coefficients from the 
nominal case. 


intended to test the feature of our estimation procedure that enables it to operate on time samples from a 
continuous nonlinear model. Fig. 6 shows the estimation results for the linear lift and drag coefficients after 
applying the algorithm separately to the nominal and damage cases. The same broad prior distribution was 
used for each case, and the algorithm processed each data set of 450 samples in about 100 ms converging 
in only 4 iterations of the main loop. These results demonstrate the algorithm’s ability to recover the 
actual coefficients and quantify uncertainty. In this case, the damage is significantly more distinguishable 
in the linear lift coefficient than in the drag coefficient. Note also that the algorithm adaptively determines 
estimation uncertainty based on the data it receives, and in this example the measured data enables a more 
accurate quantification of the damage model coefficients. 

So far our analysis has only examined the ability to discern the fault after knowing when it occurred. To 
demonstrate ability to identify the fault we turn to our evidence calculation, and instead of performing bulk 
estimates after all the data is collected, our estimation algorithm is executed each time a new observation 
appears. Furthermore, each estimate is obtained using only the 20 most current observations (i.e., 2 s of the 
most recent data). The results of this calculation are plotted in Fig. 7 for two different prior assumptions. 
The first, nominal prior case, assumes that the nominal coefficients are known to high precision ahead of 
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Figure 6: Estimation results from trajectory and accelerometer measurements for nominal and 20% 
fault models. 


time, corresponding to the case where one expects the aircraft to behave normally. The second, open prior 
case, assumes the prior is a broad distribution on the model parameters, corresponding to an openness to 
learn the behavior of the aircraft from only its most recent maneuvering performance. These two cases are 
viewed as two separate modeling assumptions with different evidence calculations. At the onset of the fault, 
the evidence immediately collapses in both cases. The calculation then remains collapsed for the strong high- 
precision prior, because the measured data contradicts the prior assumption in this case. This is in contrast 
to the open prior case, which recovers the faulted aerodynamic coefficients and evidence level within the 20 
sample window. The evidence collapse in each case is most significant at the onset of the fault, because in 
each case the new faulted data point contradicts the current state of knowledge given the past observations 
most strongly. To examine this phenomenon a little closer, the bottom plot shows a histogram of the change 
in evidence as each new measurement is added (and an old measurement is forgotten). The change caused 
by the fault is clearly several standard deviations from the natural variation in evidence for our small sample 
window, demonstrating a near instantaneous ability to detect the fault in this case. 

IV. Maneuverability Beyond Trim 

Based on the aerodynamic information recovered through the system identification process discussed in 
§111, and on the trim envelope calculated in §11. B, one is now able to obtain extended safe maneuvering 
envelopes by using the robust reachability analysis summarized in this section. 

IV. A. Interpretation of the Maneuvering Envelope 

In this context, the preferred interpretation of the safe maneuvering envelope considers reachability from the 
trim envelope. The stable and controllable trim envelope is considered an a-priori safe set. The backwards 
reachable set is defined as the set of states from where (at least one point in) the trim envelope can be 
reached. The forwards reachable set is defined as the set of states which can be reached from (at least 
one point in) the trim envelope. The safe maneuvering flight envelope is then the intersection between the 
forwards and backwards reachable sets. This interpretation is illustrated in Fig. 8. In addition to the safe 
envelope, the backwards reachable set is considered as the survivable flight envelope. After an upset due to 
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Figure 7: Fault identification from evidence. The top plot shows the running evidence calculation, 
for two separate prior assumptions, based on a 20 sample window of the most recent data. The bottom 
plot shows a histogram of the change in evidence as each new sample is added into the window (and 
the last sample is forgotten) . Clearly, the change in evidence caused by the fault is significantly greater 
than the natural variation in evidence caused by new measurements. This enables near instantaneous 
identification of the fault in this scenario. 


damage, turbulence, a wake encounter etc., it is possible to bring the aircraft back to a safe trim condition 
as long as the current flight condition is situated inside the backwards reachable set. 

The goal, then, is to perform a combined forward and backward reachability analysis from the trim 
envelope as efficiently as possible, in an effort to enable online implementations. Based on previous re- 
search, 37 level set methods are an excellent candidate. Finally, robustness to the uncertainty in the system 
identification (i.e., aerodynamic parameter estimates) is an important aspect to consider as well. 

IV. B. Optimal Control Formulation 

It has been shown in the literature that maneuvering envelope estimation through reachability can be refor- 
mulated in the optimal control framework. 9 Consider our continuous time control system: 

x = /(x, u ; c* + A) (34) 

with (x, u) G 6, as defined in §11. B above. Here we consider the case where the aerodynamic coefficient vector 
can take on any value within a probable set {c = c* + A | A G D C M d }, such as a 95% confidence ellipse. 
Define 4> (t, £, x, u (•) , A) as the state trajectory, where A characterizes possible parameter uncertainty. Given 
a set of states if C M n , and a set of possible input functions the reachability question is naturally 
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Figure 8 : Safe maneuvering envelope as intersection between forwards and backwards reachability, 
modified from source: van Oort. 12 


formulated using the relation between the set K and the state trajectories 0 of (34) over the horizon T. The 
problem of interest is the following: 

Robust reachability: Does there exist a u G W[ 0 j t] and a f G [0,T] such that the trajectory </> of the 
state x satisfies x G K, irrespective of A? 

The answer to this question is found by playing a pursuit evasion game over the horizon T > 0 with target 
set K C M n . 38 It is assumed that u is trying to bring or keep the state in the set K, whereas A is trying to 
drive it out of K. To ensure the game is well-posed, u is restricted to play non-anticipative strategies with 
respect to the parameter uncertainty A E D. 

For the types of safety problems considered here, a set of initial states has to be established such that u 
can win the game. To this end, the robust reachable set 1Z a is defined as, 


lZ A (t, K) = {x G K n | VA G D, 3u eU [t> T], [t,T], , A) e K) . 


(35) 


As shown in the literature, a useful characterization of this set is found by solving an INFMIN problem . 9 
The crux is to include the A’s as disturbances opposing the optimization over u. Consider a closed set K, 
defined as the level set of a continuous function l : M n -A M, i.e., K = {x G M n | l (x) > 0}. We then define 
the invariance set as , 11 

In v(t,K) = {xe M n | V 2 (x,t) > 0 } , (36) 


where, 


V 2 (x,t) = inf sup min l(<p(r,t,x,u(-), A)). 
u(-)eU [tjT ] A £D re[t,T\ 


(37) 


The solution for V 2 is then found by solving the associated Hamilton- Jacobi-Bellman (HJB) Partial Differ- 
ential Equation (PDE ): 9,38 


dV 2 

dt 


f dV 2 

(x, t) + min < inf sup — — 

re[t,T\ [u(-)eu [ttT] AeD ox 


(x,t)f(x,u, A) \ = 0, 


(38) 


where V 2 (x,T) = l (x) holds for backward integration and V 2 (x,t) = l (x) applies to forward integration. 
These HJB PDE’s can be solved using level set methods, for which a toolbox is available in Matlab ®. 10 
Finally, the robust reachability set is found through duality: 1l/±(t,K) = [Inv(£, K c )] c . See Ref. [39] for 
further details. 


IV. C. Examples 

The scenario to be considered is the (V, 7 ) maneuvering envelope for the previously introduced RCAM model 
with bank angle ip = 0. Both the nominal and generic damage scenarios will be evaluated. The damage 
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scenario is the one considered previously in §11. C, with a 20% decrease (increase) in lift (drag) performance, 
as well as envelope degradation caused by 50% loss of maximum thrust authority. This example builds further 
on the results for the trim envelope presented in § II. C. The results are shown in Figs. 9-11. Backward 
reachability analyses corresponding to the survivable maneuvering envelope are presented in Fig. 9. Fig. 
10 shows the forward reachability analysis results. The intersection of forwards and backwards reachability 
is defined as the safe maneuvering envelope and is shown in Fig. 11. The time horizon in this example is 
set at T = 5 s. This choice is based on the time constants of the considered relevant dynamics. The blue 
rectangular contour corresponds to the largest rectangle which can be drawn in the trim envelope of the 
nominal aircraft as depicted in Fig. 4. Similarly the red contours correspond to the largest rectangles which 
can be spanned in the trim envelopes of both damage scenarios in Fig. 4. 

Fig. 9 analyzes backward reachability from the trim envelopes over a time horizon of T = 5 s. By 
comparing the undamaged and damaged envelope boundaries, one can see the influence of the damage 
characteristics. It can be seen that for lower speed conditions, recovery with damage is only possible for 
a more restricted range of negative flight path angles. The increased drag results in a higher minimum as 
well as maximum airspeed limit, since accelerating works slower and decelerating works faster. The loss in 
lift capability means that smaller positive flight path angles are backwards reachable at lower airspeed. The 
scenario with the reduced maximum thrust limit is similar, but negative flight path angles at lower speeds 
are even more restricted due to the extra limitation in thrust. Furthermore, the reduction in backwards 
reachable positive flight path angles is more significantly influenced by the reduced trim envelope, than by 
the actual reachability calculations over 5 s. 

Comparison of backwards reachability, damaged and undamaged Comparison of backwards reachability, damaged and undamaged 




V [m/s] V [m/s] 

(a) Effect of a 20% decrease in the lift coefficients, and a 20% (b) Additional envelope degradation caused by 50% loss of 
increase in the drag coefficients maximum thrust authority 

Figure 9: Backward reachability analysis over T — 5 s, based on calculated trim envelope boundaries 
and identified aerodynamic parameters. Upper left envelope areas correspond to nominal configuration. 

Fig. 10 shows the forward reachable set calculations over the same time horizon. One may observe that 
there is less range across larger flight path angle values in the damage scenarios due to the reduction in lift 
capability. Fig. 10(a) shows that larger positive flight path angle values are still reachable at higher speeds 
with damage. This is related to the fact that the trim envelope of the damaged scenario has moved towards 
a higher airspeed range, which makes comparing both envelopes less straightforward. Forwards reachability 
for the scenario depicted in Fig. 10(b) is similar, but forwards reachability for positive flight path angles is 
further reduced due to the changed trim envelope boundaries and the extra limitation on thrust. Comparing 
Fig. 10(a) and 10(b), one finds that the reduced thrust authority also affects the reachable speed range, 
which is now more restricted at the upper and lower boundaries. 
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Comparison of forwards reachability, damaged and undamaged 



V [m/s] 


Comparison of forwards reachability, damaged and undamaged 



V [m/s) 


(a) Effect of a 20% decrease in the lift coefficients, and a 20% (b) Additional envelope degradation caused by 50% loss of 
increase in the drag coefficients maximum thrust authority 


Figure 10: Forward reachability analysis over T = 5s, based on calculated trim envelope boundaries 
and identified aerodynamic parameters. Upper left envelope areas correspond to nominal configuration. 


The safe maneuvering envelopes presented in Fig. 11 are the result of taking the intersection between the 
backward and forward reachability envelopes in Figs. 9 and 10. In Fig. 11(a), it can be seen that the safe 
maneuvering envelope shift is similar in nature to the trim envelope shift shown in Fig. 4. A higher speed 
range is needed to compensate for the loss in the lift force coefficient, and smaller positive flight path angles 
can be reached because smaller excess of net thrust is available with increased drag. In Fig. 11(b), one may 
observe that the restriction on maximum thrust leads to a further shrinking of the maneuvering envelope, 
especially for positive flight path angles. Level and climbing flight are still possible, but only in non-trim 
and as a trade-off between kinetic and potential energy. One has to sacrifice speed for non-negative flight 
path angles. In practice, however, it is the other way around: one would rather sacrifice flight path angle 
for airspeed, in order to prevent stalling. 

Extensive Monte Carlo analyses were also performed in order to verify the accuracy of the boundaries 
of the estimated maneuvering envelopes. These analyses were based on the non-simplified aircraft model, 
ignoring the assumption that the aerodynamic angles a and /3 should be small. All these Monte Carlo 
analyses have confirmed that the results provided here are accurate and that the simplifications hold for the 
current ranges of the aerodynamic angles, namely a 6 [0°;14.5°] (no stall) and (3 E [—5°; +5°]. This is an 
important conclusion which makes a relevant online safe maneuvering envelope estimation tool much more 
feasible. This result is covered in [39], along with more examples and discussion. 

V. Major Application Areas 

V.A. Flight Management Systems 

The Flight Management System (FMS) of commercial aircraft is a common tool to aid crew in piloting 
aircraft. The FMS is comprised of many functions, including navigation, flight planning, trajectory genera- 
tion, performance computations, and guidance. With the exception of navigation, all these capabilities are 
either directly or indirectly dependent on aircraft performance characteristics; and when these characteris- 
tics change due to vehicle impairment, the ability of the aircraft to satisfy the constraints of the existing 
flight plan or to follow the recommended trajectory, guidance commands, or performance profiles is compro- 
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(a) Effect of 20% decrease in lift coefficient and 20% increase 
in drag coefficient 
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(b) Additional envelope degradation caused by 50% loss of 
maximum thrust authority 


Figure 11: Calculation of safe maneuvering envelope sets through forward and backward reachabil- 
ity analysis over T = 5s, based on calculated trim envelope boundaries and identified aerodynamic 
parameters. Upper left envelope areas correspond to nominal configuration. 


mised. Therefore, the trim and maneuverability envelope calculations discussed in this paper enable envelope 
protection enhancements for next generation flight management systems, like the one discussed in [40]. 

V.B. Pilot Cockpit Displays 

Pilot cockpit displays are instrumental for pilot situational awareness of aircraft state and performance 
limitations. Performance limitations such as V-speeds or maximum and minimum speeds for the nominal 
configuration are currently determined from performance databases that will become obsolete in the event 
of failure, damage, or impairment conditions. Even under nominal conditions, the relationship between 
maneuverable airspeed and vertical speed changes in non-obvious ways with altitude, air density, aircraft 
configuration (e.g., flap settings, stabilizer position), and bank angle. A recent survey of airline transport 
rated pilots demonstrated a preference for the display of maneuvering limitations in bank angle, vertical 
velocity, and aircraft speed, 41 all of which are readily provided by the envelope characterization methods 
discussed here. 


VI. Conclusion 

There are three distinguishing features of the maneuverability characterization methods presented in this 
article. 

First, instead of linearizing the full aircraft dynamics model, we derived a useful nonlinear model that 
allowed us to avoid local linear approximations all together. This enabled rapid numerical computation 
of the curved (non-convex) trim envelope, which is more representative of the global aircraft performance 
capability than would be obtained through the use of linearized models. 

Second, our integrated system identification approach leverages the same dynamics model to estimate the 
underlying aerodynamic parameters from currently available flight data, while also establishing confidence 
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regions biased towards greater uncertainty when insufficient input excitation is available. This is impor- 
tant because the penalty for over- confidence in-flight can be severe, while the cost of under- confidence is 
often only over-preparedness and caution. Furthermore, the system identification process seamlessly pro- 
vides methods for learning the aircraft dynamics passively in the background, detecting impairment, and 
also selecting between multiple nonlinear dynamics models for damaged aircraft through the calculation of 
Bayesian evidence. 

Third, through the safe maneuvering envelope calculation, we effectively maximized the options pilots 
or automation systems have to recover the aircraft in damage scenarios, while remaining robust to the 
uncertainty in the system identification process. 

The net result is an increased flexibility for developing advanced aircraft diagnostics that provide the 
bottom line maneuverability of the aircraft as an output, and this is expected to have important applications 
to flight planning, trajectory generation, guidance algorithms, and pilot displays. 
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